\documentclass[preprint,amsmath,pre,aps,superscriptaddress,showpacs]{revtex4}
%\documentclass[a4paper,11pt]{repscr_pc}
\usepackage{texdraw}
\usepackage{graphicx}
\usepackage{epsf}
\usepackage{subfigure}




\usepackage{amsmath}
\usepackage{float}


%\setlength {\textwidth}{16cm}
%\setlength {\topmargin}{0cm}
%\setlength {\textheight}{23cm}
%\oddsidemargin 0mm
\input{macrosLBGK}

\begin{document}
\input txdtools

\title{\bf Lattice-Boltzmann simulation of single phase and immiscible fluid flow in Bentheimer sandstone}

\author{E.\ S.\ Boek\footnote{E-mail: e.boek@imperial.ac.uk}}%, esb30@cam.ac.uk}}
\affiliation{Department of Chemical Engineering, Imperial College London, South Kensington Campus, London SW7 2AZ, United Kingdom}
\affiliation{Schlumberger Cambridge Research, High Cross,
           Madingley Road, Cambridge CB3 0EL, United Kingdom}
%\affiliation{Schlumberger Cambridge Research, High Cross,
%           Madingley Road, Cambridge CB3 0EL, United Kingdom}
\author{P.\ V.\ Coveney}
\affiliation{Centre for Computational Science, Christopher Ingold
Laboratories, \\ University College London, 20 Gordon Street, London
WC1H 0AJ, UK}
\author{J.\ Harting}
\affiliation{Department of Applied Physics, TU Eindhoven,
Den Dolech 2, 5600 MB Eindhoven, The Netherlands}
\author{M.\ Venturoli \footnote{current address: Courant Institute for Mathematical
Sciences, NYU, 251 Mercer Street, New York 10012 N.Y., U.S.A.}}
%E-mail: mventuro@cims.nyu.edu}}
\affiliation{Schlumberger Cambridge Research, High Cross,
           Madingley Road, Cambridge CB3 0EL, United Kingdom}
%\affiliation{Courant Institute for Mathematical
%Sciences, NYU, 251 Mercer Street, New York 10012 N.Y., U.S.A.}
\affiliation{Centre for Computational Science, Christopher Ingold
Laboratories, \\ University College London, 20 Gordon Street, London
WC1H 0AJ, UK}

\date{\today}

%______Abstract
\begin{abstract}
        We use lattice-Boltzmann (LB) simulations with multi-relaxation time (MRT) collision term to calculate the flow in a
complex three-dimensional porous medium, in this case Bentheimer
sandstone. The data sets used in this study were obtained using
X-ray microtomography (XMT) on a cylindrical core of Bentheimer
sandstone. First, we consider single phase flow. We systematically
study the effect of system size and periodic boundary conditions and
validate Darcy's law from the linear dependence of the flux on the
body force exterted. Then, we find that, as the system size of the
computational sub-sample increases, the values of the permeability
measurements as a function of porosity tend to concentrate in a
narrower region of the porosity. We observe that for different
boundary conditions the permeabilities of the larger sub-samples are
less dependent on the particular boundary conditions used. Then we
calculate displacement propagators using the LB flow field and
observe dispersive behaviour, in agreement with NMR experiments.
Finally, we consider immiscible flow and calculate relative
permeabilities for different forcing conditions.
\end{abstract}

\pacs{}

\maketitle


%______Introduction
\section{Introduction}
For many scientific and industrial applications, it is important to
know the permeability of a particular porous medium, such as
sandstone. One important application is in the oilfield industry,
where the production of hydrocarbons is often limited by the
permeability of the reservoir rock. Recently, lattice-Boltzmann (LB)
simulations have become available as a method to calculate the
permeability of complex three-dimensional geometries. These
simulations also give access to the full flow field at the pore
scale. For single phase flow, the LB flow field generated can then be
used to calculate hydrodynamic dispersion and displacement propagators
in the porous medium. For the case of multi-phase flow, it is very
cumbersome to reliably determine the relative permeability of a rock
core using experimental methods. For this reason, the LB method has
become popular as a realistic alternative to calculate the relative
permeability of complex three-dimensional geometries. A number of
studies, using the LB methodology, has been reported in the
literature.
Ferr\`eol and Rothman were among the first to use LB
simulations to study flow through Fontainebleau
sandstone~\cite{bib:ferreol-rothman}.
Auzerais {\em et al.}~\cite{bib:auzeraisGRL96} used LB
simulations to calculate the permeability of Fontainebleau
sandstone. They used X-ray Micro-Tomography (XMT) imaging and provided
a direct comparison with laboratory measurements.
Kang {\em et al.}~\cite{Kang-2002} developed a unified LB
method for flow in multiscale porous media. This model not only can
simulate flow in porous systems of various length scales but also can
simulate flow in porous systems where multiple length scales
coexist. Simulations performed on a fractured porous system show that
this method gives very good overall permeability values for
the whole fractured system.
Martys and Chen~\cite{bib:martys-chen} described a LB method to
study multi-component flow in complex three-dimensional porous media.
Sukop {\em et al.}~\cite{bib:sukop} studied multi-phase fluids in
porous media which were compared with X-ray micro-tomography.
Keehm {\em et al.}~\cite{bib:keehm-2004} studied the
permeability of rock samples using 3D reconstruction and LB flow
simulation.
Although the list above is far from exhaustive, the papers mentioned
provide a good introduction to the subject. In this paper, we report
on the calculation of the flow in a realistic complex
three-dimensional (3d) sample of Bentheimer sandstone.
In section~\ref{Sec:LatticeBoltzmann} we briefly introduce the
lattice Boltzmann method. The porous medium and the quantity computed are
described in section~\ref{sec:material}.
Results on the dependence of the permeability measurement
on the system size of the computational sub-sample are reported
in section~\ref{sec:results}.
%With increasing
%sub-sample size, we find that the values of the permeability as a
%function of porosity tend to concentrate in a narrower region of the
%porosity. %already said in abstract and conclusions
There we also investigate the effect of boundary conditions.
We observe that for different boundary conditions, the permeabilities of
the larger sub-samples are less dependent on the particular boundary
conditions used than those of the smaller sub-samples.
In section~\ref{propagators}
we calculate displacement propagator distributions using the LB flow
field. Finally, in section~\ref{relperm},
we calculate relative permeability curves for the case
of immiscible flow.
The current paper follows up on a recent study, where we reported preliminary results
of a LB study of fluid flow in two- and three-dimensional
porous media with realistic rock geometries~\cite{bib:boek-venturoli-2009}.


%______Lattice Boltzmann theory
\section{Lattice-Boltzmann method}
\label{Sec:LatticeBoltzmann}

During the last decade, many authors have shown that the
lattice-Boltzmann algorithm is a powerful method for simulating fluid
dynamics.  This success is due to its simplicity and to facile
computational implementations
\cite{bib:chin-harting-jha,bib:love-nekovee-coveney-chin-gonzalez-martin,bib:nekovee-chin-gonzalez-coveney,bib:succi}.
Instead of tracking individual atoms or molecules, the
lattice-Boltzmann method describes the dynamics of the single-particle
distribution function of mesoscopic fluid packets.

In a continuum description, the single-particle distribution function
$\fone$ represents the density of fluid particles with position
$\vect{r}$ and velocity $\vect{v}$ at time $t$, such that the density
and velocity of the macroscopically observable fluid are given by
$\rho(\vect{r},t) = \int \fone {\mathrm d}\vect{v} $ and
$\vect{u}(\vect{r},t) = \int \fone \vect{v} {\mathrm d} \vect{v}$
respectively. In the non-interacting, long mean free path limit, with
no externally applied forces, the evolution of this function is
described by the Boltzmann equation
\begin{equation}
  \label{eq:boltzmann}
  \left( \dt + \vect{v} \cdot \vect{\nabla} \right) f_1
  = \Omega[f_1].
\end{equation}
While the left hand side describes changes in the distribution
function due to free particle motion, the right hand side models
pairwise collisions. The collision operator $\Omega$ is an integral
expression that is often simplified to the linear
Bhatnagar-Gross-Krook (BGK)~\cite{bib:bgk} form
\begin{equation}
  \label{eq:bgk}
  \Omega[f] \simeq - \frac 1 \tau \left[ f -
    f^{\mathrm{(eq)}} \right].
\end{equation}
This collision operator describes the relaxation, at a
rate controlled by a characteristic time $\tau$, towards the local
Maxwell-Boltzmann equilibrium distribution $f^{\mathrm{(eq)}}$. It can
be shown that distributions governed by the simple Boltzmann-BGK
equation conserve mass, momentum, and energy \cite{bib:succi}. They
obey a non-equilibrium form of the second law of
thermodynamics~\cite{bib:liboff} and the Navier-Stokes equations for
macroscopic fluid flow are obeyed on coarse length and time scales
~\cite{bib:chapman-cowling,bib:liboff}.

By discretising the single-particle distribution function in space and
time, one obtains the lattice-Boltzmann regular formulation, where the
positions $\vect{r}$ on which $\fone$ is defined are restricted to
points $\vect{r}_i$ on a Bravais lattice. The velocities $\vect{v}$
are restricted to a set $\ci$ joining points on the lattice and the
density of particles at lattice site $\vect{r}$ travelling with
velocity $\ci$, at timestep $t$ is given by $f_i(\vect{r},t) =
f(\vect{r},\ci,t)$.  The fluid density is given by
%
\begin{equation}
  \label{eq:lbe-density}
  \rho(\vect{r}) = \sum_i f_i(\vect{r}),
\end{equation}
and the velocity  by
\begin{equation}
  \label{eq:lbe-velocity}
  \vect{u}(\vect{r}) = \sum_i f_i(\vect{r}) \ci.
\end{equation}
%
The discretised Boltzmann distribution function can be evolved in a two-step
procedure. In the collision step, the densities at each lattice site are
redistributed across the velocity vectors; this process corresponds
to the action of the collision operator. In the advection step,
values of the post-collisional distribution function are propagated
to adjacent lattice sites.

By combining the two steps, one obtains the lattice-Boltzmann equation
(LBE)
%
\begin{equation}
  \label{eq:lbgk}
  f_i(\vect{r}+\ci \Delta t ,t+\Delta t) - f_i(\vect{r},t)
  = \Omega[f] \\
  = - \frac 1 \tau \left[ f_i(\vect{r},t)
    - N_i\left( \rho, \vect{u} \right)\right],
\end{equation}
%
where $N_i = N_i\left(\rho(\vect{r}),\vect{u}(\vect{r})\right)$ is a
polynomial function of the local density and velocity, and can be
found by discretizing the Maxwell-Boltzmann equilibrium distribution.

To model multicomponent interacting fluids, our lattice-Boltzmann
implementation uses the Shan-Chen approach~\cite{bib:shan-chen}, which
incorporates an explicit forcing term in the collision operator.  Shan
and Chen extended the single-particle distribution function $f_i$ to
the form $f_i^\sigma$, where $\sigma$ denotes different fluid components.
The density and momentum of component
$\sigma$ are given by $\rhosig = \sumi \fis$ and $\rhosig
\vect{u}^\sigma = \sumi \fis \ci$ respectively. The fluid viscosity
$\nu^\sigma$ is proportional to $(\tau^\sigma-1/2)$. This results in a
lattice BGK equation for each component $\sigma$ of
the form
\begin{equation}
  \label{eq:lbgk-sc}
  \fis(\vect{r}+\ci \Delta t ,t+\Delta t) - \fis(\vect{r},t) =
  - \frac 1 \tausig
  \left[
    \fis - N_i(\rhosig, \vect{v}^\sigma)
  \right].
\end{equation}
The velocity $\vect{v}^\sigma$ is found by calculating a weighted average
velocity
\begin{equation}
  \vect{u}' = \left( \sumsig \frac \rhosig \tausig \vect{u}^\sigma
  \right)
  / \left( \sumsig \frac \rhosig \tausig \right),
\end{equation}
and then adding a term $\vect{F}^\sigma$ to account for the additional forces
\begin{equation}
  \vect{v}^\sigma = \vect{u}' + \frac \tausig \rhosig \vect{F}^\sigma .
\end{equation}
To produce nearest-neighbour interactions between fluid
components, this term assumes the form
\begin{equation}
  \label{eq:colour-colour}
  \vect{F}^\sigma =
  - \psisig ( \vect{r} )
  \sumsigb g_{\sigma \bar{\sigma}}
  \sumi \psisigb \left( \vect{r} + \ci \right) \ci,
\end{equation}
where $\psisig ( \vect{r} ) = \psisig ( \rhosig ( \vect{r}))$ is an
effective charge for component $\sigma$.  The easiest choice for
$\psisig ( \vect{r} )$ is to set it equal to the density of the fluid,
i.e. $\psisig ( \vect{r} ) = \rhosig ( \vect{r})$.  This is the
functional form we have chosen for the simulations presented in this
work.  $g_{\sigma \bar{\sigma}}$ is a coupling constant controlling
the strength of the interaction between two components $\sigma$ and
$\bar{\sigma}$. If $g_{\sigma \bar{\sigma}}$ is set to zero for
$\sigma = \bar{\sigma}$, and to a positive value for $\sigma \neq
\bar{\sigma}$ then, in the interfacial region between bulk domains of
each component, particles experience a force in the direction away
from the interface, producing immiscibility. For two-component
systems, we use the notation $g_{\sigma \bar{\sigma}} =
g_{\bar{\sigma}\sigma} = g_{br}$, where $b$ and $r$ denote the
``colour'' of the particles (blue/red, respectively).  External forces
are added in a similar manner. For example, in order to produce a
gravitational force acting in the $z$-direction, the force term
$\vect{F}^\sigma$ can take the form $g \rho^\sigma \hat{\vect{z}}$.

A convenient way to characterize binary fluid mixtures is the order
parameter or colour field
\begin{equation}
  \phi(\vect{r}) =
  \rho^r(\vect{r})-\rho^b(\vect{r}).
  \label{Eq:Colour}
\end{equation}
The order parameter is positive in areas of high concentration of
`red' fluid and negative in areas of `blue' dominance; the isosurface
$\phi(\vect{r}) = 0$ denotes the interface between the two fluid
constituents.

This model has also been extended to handle amphiphiles, which are
treated as point dipoles with different interaction
strengths on each end and a rotational degree of
freedom~\cite{bib:chen-boghosian-coveney}. For the simulations we
use the code LB3D~\cite{lb3d}, which can handle binary and ternary
fluid mixtures with or without amphiphiles. Since here we only
discuss simulations of binary fluids, we refer the reader to
references
\cite{bib:chen-boghosian-coveney,bib:nekovee-chin-gonzalez-coveney,bib:love-nekovee-coveney-chin-gonzalez-martin}
for a more comprehensive description of the amphiphilic case.
%\fixme{
%A well-known drawback of the lattice-Boltzmann method is its numerical
%instability which occurs if too high forcing rates or too high interparticle
%interactions are applied.
%}


\subsection{Multi-Relaxation-Time (MRT) scheme for the lattice Boltzmann method}
MRT scheme allows independent adjustment of bulk

%\footnote{Bulk viscosity, also called volume viscosity is important only for the simulations where fluid compressibility is essential. It appears in the compressible Navier-Stokes equation:  
%\begin{equation}
%\rho(\frac{\partial v}{\partial t}+v \cdot \nabla v)=-\nabla \rho +\mu\nabla^2 v+f+\mu^v\nabla(\nabla \cdot v)
%\end{equation}
%}
%where $\mu^v$ is the bulk viscosity. In the incompressible Navier-Stokes equation, this term disappears because the divergence of the velocity of an incompressible fluid, $\nabla \cdot v$, equals 0 } 
and shear viscosities which significantly improves the numerical stability for a low viscosity fluid. In the single relaxation time LB model, the collision term is relaxed by a single parameter $\tau$, while it could be relaxed instead by matrixs:

\begin{equation}
        f_i(x+c_i \Delta t)-f_i(x,t)=M_{-1}S((Mf)-m^{eq}
\end{equation}

The transformation matrix $M$ is comprising of constants. The value of components can be found in \cite{lallemand2000}:



The macroscopic quantities are calculated in the same way as in the LBGK model. Instead of considering distribution functions, MRT employs several moments corresponding to macroscopic quantities and their flux. These quantities can be relaxed with different time scales. 


The vector $m^{eq}$ is the equilibrium moments and are given by \cite{tolke_2006}:
\begin{eqnarray}
        &m^{eq}_0=\rho,&\\
        &m^{eq}_1= e_{eq}=\rho_0(u_x^2+u^2_y+u^2_z),&\\
        &m^{eq}_3=\rho_0 u_x,&\\
        &m^{eq}_5=\rho_0 u_y,&\\
        &m^{eq}_7=\rho_0 u_z,&\\
        &m^{eq}_9=rho_0(2u^2_x-u^2_y-u^2_z),&\\
        &m^{eq}_11=\rho_0(u^2_y-u^2_z),&\\
        &m^{eq}_13=rho_0u_xu_y,&\\
        &m^{eq}_14=rho_0u_yu_z,&\\
        &m^{eq}_15=rho_0u_xu_z,&\\
        &m^{eq}_2=m^{eq}_4=m^{eq}_6=m^{eq}_8=m^{eq}_16=m^{eq}_17=m^{eq}_18=0,&
        \end{eqnarray}

The matrix "S" is a diagonal matrix with relaxation parameters:

\begin{equation}
        S=(0,s_e,s_{e^2},0,s_q,0,s_q,0,s_q,s_{\nu},s_{\pi},s_{\nu},s_{\pi},s_{\nu},s_{\nu},s_{\nu},s_{m},s_m,s_m)
\end{equation}

The kinematic viscosity is given by:
\begin{equation}
\nu=c^2_s(\frac{1}{s_{\nu}}-\frac{1}{2})\Delta_t 
\label{nu_cal}
\end{equation}

where $c_s$ is the local sound speed that equals $1/\sqrt{3}$ in the single phase LB model, the other parameters can be chosen freely by users. In this paper we use a choice of these parameters  suggested by \cite{Pan_2006} which eliminate the dependence of permeability calculation on $s_{\nu}$:

\begin{equation}
        s_e=s_{e^2}=s_{\pi}=s_{\nu}, s_q=s_m=8\frac{2-s_{\nu}}{8-s_{\nu}}
        \end{equation}




\subsection{Boundary conditions at the solid nodes: a note of warning}
\label{sec:BC}

To study flow in porous media at the pore scale using LB simulations,
it is very important to use appropriate boundary conditions for
fluid-rock interactions.  Wall boundary conditions in the LB method
were originally taken from the lattice gas (LG) method.  For example,
the bounce-back scheme is used at the walls to obtain no-slip velocity
conditions.  In the bounce-back scheme, when a particle distribution
streams to a wall node, it scatters back to the node it came from.
%
Because for a node near a boundary some of its neighboring nodes may be
located outside the flow domain, the distribution functions from these
non-fluid nodes are therefore unknown after each streaming process.
The bounce-back scheme is a simple way to fix these unknown
distributions on the wall nodes. It was found that the bounce-back
condition is only first-order in numerical accuracy at the boundaries
\cite{cor911,zie931,gin941}. However, if the location of the wall is
placed half-way between the wall node and the first adjacent fluid
node, it can be shown that the method is second order.
%
These results were proven both by numerical experiments and by the
analytical solution of the LBGK equation for a Poiseuille flow in 2d
\cite{zou951,bib:he-zou-luo-dembo}.  However, even with the half-way
bounce-back boundary condition, the Poiseuille profile shows
deviations from the analytical solution, i.e., it is still a
parabola but with a shift at the boundary.  The bounce-back boundary
condition mimics boundaries that move with a speed that depends on
the relaxation time $\tau$.  In the general case, the precise
location of the wall at the zero slip velocity lies somewhere
between the wall node and the neighboring fluid node and can be
determined by the eigenvalues of the linearised collision operator
\cite{gin941}.

%PVC:
It should be noted that, in this paper, we will be using a single relaxation time $\tau$. The great advantage is that this scheme can be easily implemented. However, this does yield a quasi-incompressible fluid whose behaviour coincides with solutions of the incompressible Navier-Stokes equations in the limit of low Mach number. We will consider the effects of the single relaxation time in the results presented in the following sections.

%Please note that in this paper, we will be using a single relaxation time $\tau$.
%The great advantage is that this scheme can be easily implemented.
%However, it comes also with a disadvantage, in the sense that it yields a compressible fluid.
%We will analyse the effects of the single relaxation time in the results presented in the following sections.


\subsubsection{A brief review of some second-order boundary conditions}
A solution to these numerical errors are the so-called second order
boundaries, for which the desired location of the non-slip
boundaries is determined by extrapolating the distribution function
from the last fluid points.

Inamuro et al. \cite{bib:inamuro} proposed a method in which the
unknown distribution functions are assumed to be an equilibrium
distribution function with a counter slip velocity, which is
determined so that the fluid velocity at the wall is equal to the
wall velocity. The limitation of this method is that it can handle
only flat boundaries and has problems with nodes at the corners.

Noble et al. \cite{nob951} introduced a method to calculate the
particle distribution at the boundaries from the velocity boundary
conditions and the particle distributions of the neighbouring fluid
nodes (utilising conservation of mass and momentum). They
seek an answer to the question: knowing the velocity of the fluid at
the boundaries, can an appropriate particle distribution be found
such that imposing this distribution provides the prescribed
velocity condition? This hydrodynamic approach seeks to maintain a
specified velocity profile on the boundaries. Moreover, the density
at the boundary becomes a computed quantity. Since the density is
related to the pressure through the isothermal equation of state,
the algorithm they propose supplies the pressure boundary condition.
However, also this method has the drawback that it can be applied only to
simple geometries, such as flat boundaries.

Chen et al. \cite{che961} enforced the condition of correct flux of
momentum near the wall as an extra condition. In this approach, the
lattice-Boltzmann scheme is regarded as a special finite difference of
the kinetic equation for $f_i$ in the Boltzmann equation for the
discrete velocity distribution.
%\begin{equation}
%  \frac{\partial f_i}{\partial t} + {\bf e}_i \cdot {\bf \nabla}f_i = \Omega_i.
%\end{equation}
These authors proposed boundary conditions as an extrapolation scheme,
to allow the boundary nodes to evolve following the lattice Boltzmann
relaxation dynamics. This is achieved by assuming that there is an
additional layer of nodes inside the wall beyond the boundaries. The
distribution functions of these additional nodes are calculated before
the streaming operation using a second order extrapolation based on
the distributions on the boundary layer and the neighbouring layer
inside the fluid. After this extrapolation, the streaming step is
carried out for all nodes. The collision step for the boundary (wall)
nodes is done by enforcing the velocity (or pressure) boundary
conditions for the equilibrium distribution functions.

Verberg and Ladd \cite{bib:verberg-ladd-01} proposed a continuous
bounce-back boundary condition to incorporate solid-fluid boundary
conditions on length scales smaller than the grid spacing. They
assumed that the population densities are uniformly distributed
throughout the volume of the Wigner-Seitz cell surrounding each
node. For each node, a real numbered parameter is then introduced,
which represents the fluid volume associated with that node.  The
distributions that stream to a solid node are reflected with rules
that depend on this fluid fraction parameter.

Martys et al. \cite{bib:martys-chen} adopted a boundary condition in
which the distribution of particles going out from the wall are set
equal to their complementary incoming velocity distributions. This
sets the normal velocity to zero. The remaining directions (in the
plane of the wall) then have their distribution functions set to the
average of the incoming directions, thus setting the tangential
velocity to zero.  During the collision step, the collision operator
is applied to the boundary nodes as well as to the fluid nodes, and
the resulting distribution functions are then propagated normally in
the next cycle \cite{zie931}.  This scheme might be a relatively
simple solution to implement for lattice Boltzmann simulations of flow
in porous media.

For practical simulations, however, the bounce-back condition is very
attractive, because it is a simple and computationally efficient
method for imposing no-flow conditions on irregularly shaped walls,
such as in porous medium.  However, since bounce-back conditions do not
give an accurate velocity, errors can be introduced in the results,
such as a permeability that depends on the fluid viscosity
\cite{Venturoli-Boek-PhysicaA}.

Another source of error in the numerical quality of LB simulations,
arises from the finite-size effects due to an insufficient number of
lattice points compared to the mean free path of the fluid particles.
These Knudsen-like effects depend on both the relaxation parameter
(which controls the mean free path) and on the lattice resolution.  The
accuracy of the results increases with increasing grid resolution, and
decreasing relaxation parameter (viscosity).

A linear dependence of the permeability on the viscosity was found in
simulations of flow in porous media that implement the bounce-back
boundary conditions
\cite{bib:ferreol-rothman,bib:koponen-kandhai-hellen-etal,bib:manwart-aaltosalmi-koponen-hilfer-timonen}.
This viscosity dependence, however, becomes weak with increasing grid
resolution.

The solutions usually adopted to solve these problems for LB studies
of flow in porous media are to use the bounce-back boundary condition,
for the sake of simplicity, and to choose a viscosity low enough to be
able to adopt a not too high spatial resolution, since the number of
lattice points increases as the third power of the resolution (for 3d
simulations).
Also, the resolution can be chosen depending on the
porosity. For large porosity, the Knudsen effects are less relevant
and the resolution can be lower than in the case of low porosity, where
the fraction of narrow channels is large
\cite{bib:koponen-kandhai-hellen-etal}.
In this study we used a value of $\nu=1/6$ ($\tau=1$), which gives accurate
results at the given resolution of the rock sample from XMT imaging.



%______Rock description and calculation of permeability
\section{Materials and methods}
\label{sec:material}
\subsection{Bentheimer sandstone}

The data sets used in this study were extracted from the image of a
cylindrical core of Bentheimer sandstone obtained by X-ray
microtomography (XMT) at the European Synchrotron Research facility in
Grenoble.  The XMT image consists of voxels arranged in a three
dimensional lattice.  The full data set consists of approximately
800x800x630 voxels, at a resolution of $\delta x$=4.9$\mu m$.  A 3d
volume rendering of the data set is shown in figure
\ref{f:BenthFullSample}.
%FIXME here add a few lines about Benthiemer: i.e. rather homogeneous, average pore/grain size, average
% permeablility, why it is important in oilfield, etc.
%----------------------------------------------------------------------
\begin{figure}[htb]
\begin{centering}
\includegraphics[width=0.44\textwidth]{./Fig/BenthFullSample.eps}
\caption{
  3d volume rendering of the XMT
  image of the Bentheimer sandstone used in this study.}
\label{f:BenthFullSample}
\end{centering}
\end{figure}
%----------------------------------------------------------------------
In the image data, each voxel corresponds to a byte-value, which
stores a greyscale value [0-255] that represents the attenuation
%reflectivity
of the rock at that point in space. This value is
proportional to the density of the material: the higher the value,
the denser is the system. Typically, the pore space has low
%reflectivity
attenuation values. The permeability and porosity of this particular
sample of rock have not been measured. However, a permeability of 2
Darcy and a porosity of 22$\%$ are typical values for Bentheimer
sandstones. % Add reference?
Using this estimated rock porosity,
it is possible to determine a threshold in the greyscale, to
discriminate between pores and rock, and byte-gate the XMT image to
generate a binary (0 pore site, 1 rock site) representation of the
rock and the pore space, as illustrated in figure \ref{f:bitmap}.
This representation can then be directly loaded into the
lattice-Boltzmann code, and used as the rock matrix (solid boundary)
for the simulations. We also define the void fraction $\phi$ in a
rock sample as the number of ``fluid'' lattice points ({\em i.e.}
the lattice points which are not on a solid obstacle) divided by the
total number of lattice points. The void fraction is related to the
porosity of the rock, which can be characterised experimentally (e.g. by
mercury porosimetry). %although a better characterisation of the
%pore space would be needed in order to define the real rock
%porosity.
%----------------------------------------------------------------------
\begin{figure}[htb]
\begin{centering}
\subfigure[greyscale]{\includegraphics[width=0.34\textwidth]{./Fig/berea.eps}}
\subfigure[bitmap]{\includegraphics[width=0.3\textwidth]{./Fig/threshold90.eps}}
\caption{
  2d slice of a subsample (128$^3$) of the
  Bentheimer sandstone used in this study.
  Figure (a) is the greyscale image from XMT data,
  and figure (b) is its bit-map representation using the value
  90 for the threshold.
}
\label{f:bitmap}
\end{centering}
\end{figure}
%----------------------------------------------------------------------

\subsection{Conversion from lattice units to physical units}
\label{s:units}
%%%
Length and time can be converted from physical units into lattice
units in the following way. The unit of length is given by the
lattice resolution $\delta x$. The unit of time, $\delta t$, can be
derived from the kinematic viscosity. We define $\nu$ to be the
kinematic viscosity in lattice units and $\nu_{\rm mKs}$ the
kinematic viscosity in mKs (SI) units.  Considering that the
kinematic viscosity has the dimensions of length squared over time,
we have
\begin{equation}
  \nu_{\rm mKs} = \nu \,  \delta x^2/\delta t.
  \label{eq:numKs}
\end{equation}
%then
%\begin{equation}
%\delta t = \nu/\nu_{\rm mKs} * \delta x^2
%\label{eq:numKs2}
%\end{equation}
The  kinematic viscosity in lattice units can be calculated using equation \ref{nu_cal}

For example, if we take the XMT resolution for the Bentheimer data
as the lattice spacing, {\em i.e.} $\delta x$=4.9$\mu m$, and the
kinematic viscosity of water, $\nu_{\rm mKs}$=10$^{-6}$
m$^2$s$^{-1}$, using equation~(\ref{eq:numKs}) with $\nu=1/6$
($s_{\nu}=1$), we have $\delta t$=4 $\mu s$.


\subsection{Calculation of rock permeabilty}
The permeability of a porous medium can be calculated from the
empirical Darcy law. This well known relation states that the flow
rate, $J$, is proportional to the force driving the fluid, the
coefficient of proportionality being the permeability of the medium,
$K$, divided by the dynamic viscosity of the fluid $\mu$. Darcy's
law can be written as
\begin{equation}
  \label{eq:darcy}
  J = -\frac{K}{\mu}(\nabla P - \rho g),
\end{equation}
where $J$ is the flow rate per unit area of cross section (flux), K is
the permeability, $\nabla P$ is the pressure drop between inlet and
outlet, $\rho$ is the fluid density, $g$ is a body force (for example
gravity), and $\mu$ is the dynamic viscosity of the fluid (related to the
kinematic viscosity by $\mu=\nu \rho$).  By measuring (or
calculating) the flux for different pressure drops (or body force
values), and using equation~(\ref{eq:darcy}), the permeability $K$ can
be derived.  The permeability has dimensions of an area, and it is
measured in units of Darcy.  Converted to SI, one Darcy becomes about
$9.869 \cdot 10^{-13}$ $m^2$, and the calculated permeability can be
converted from lattice units [l.u.] to Darcys using the following
expression
\begin{equation}
  K [{\rm Darcy}] = K [{\rm l.u.}] \frac{\delta x ^2}{0.9869}
\end{equation}
where $\delta x$ is the lattice spacing in microns.

\subsubsection{Computational details}
To estimate the single phase permeability of the sample, we impose a
flow in the positive $z$ direction of the rock.
The flow is driven only by a body force $g$
(i.e. no pressure drop is explicitly present).
This choice is due to the fact that it is difficult to
explicitly introduce a pressure difference in lattice-Boltzmann
simulations. A pressure difference could in principle be simulated by
introducing a density gradient in the fluid, but having density
differences in an incompressible fluid is not consistent. However, we
would like to point out that this approach is used by several authors,
and that it gives results in good agreement to the ones obtained by using a
body force. A correspondence between the body force, $g$, and the
pressure drop, $\nabla P$, can be defined using the following
equation,
\begin{equation}
  \nabla P=\frac{(P_{i}-P_{o})}{L}= g \rho
  \label{eq:gdeltaP}
\end{equation}
where $P_{i}$ and $P_{o}$ are the pressures at the inlet and outlet
respectively, $L$ is the distance between inlet and outlet, and $\rho$
is the fluid density.  Equation~(\ref{eq:gdeltaP}) can be used to
compare the simulations with experiments, in which usually a pressure
drop is used to drive the fluid flow.  Using the unit of length and
time, and the value of the kinematic viscosity, as described in
subsection \ref{s:units}, we have
\begin{equation}
  \nabla P = g \frac{\delta x}{\delta t^2} \, \rho \; [{\rm Pa/m}]
\end{equation}
where $\rho$ can be taken as the density of water at room temperature
(10$^3$ Kg/m$^3$), and $g$ is in lattice units.

Each simulation is run until the steady state is achieved ({\em i.e.}
until the average flux is constant in time). Then the flux across each
slice in the $xy$-plane (perpendicular to the direction of the flow)
is calculated according to
\begin{equation}
  J(z)=\frac{1}{n_x\,n_y}\,\sum_{x,y=1}^{n_x,n_y}
  v(x,y,z), \hspace{1cm} z=1,...,n_z
  \label{eq:fluxz}
\end{equation}
where $x$, $y$ and $z$ denote the coordinates of a
lattice site, and $n_x$, $n_y$ and $n_z$ are the system sizes (in
lattice points).

For spatially periodic media, the pressure can be decomposed into a
spatially periodic component $P_{per} (x)$ and a linear component $\nabla P_
x$, where $\nabla P$ is the average pressure gradient and $x$ is the
position vector. Please note that a medium which is laterally sealed
is spatially periodic along the two transversal axes.
For more details, see ~\cite{AJQ90}.
%P.M. Adler, I.C.G. Jacquin and J.A. Quiblier,
%Int. J. Multiphase Flow Vol. 16, No. 4, pp. 691-712, 1990



%%%%%%% Single phase flow %%%%%%%
\section{Results}
\label{sec:results}
\subsection{Validation of Darcy's law}
Because of mass conservation and fluid incompressibility, the flux in
each longitudinal slice parallel to the main flow direction ($z$-axis)
should be constant.  However, small compressibility effects are
observed in our lattice Boltzmann simulations, resulting in a
variation of the values of $J(z)$ with $z$, as shown in figure
\ref{f:fluxz}.
%----------------------------------
\begin{figure}[thb]
\centering
\includegraphics[width=0.6\textwidth]{./Fig/fluxz}
\caption{
  Flux, $J(z)$, across each cross sectional area
  perpendicular to the forcing---and main flow---direction ($z$-axis)
  with body force $g$=0.0001.
  %({\bf Maddalena, is this correct?}:
  %I think it is but I can't verify since I don't have the original data anymore
  %anyway, since this is the value we use in the following it seems reasonable that I have shown
  % results for that value.
  The flux is expressed in lattice units. }
\label{f:fluxz}
\end{figure}
%----------------------------------
The compressibility effects depend on the fluid velocity, and hence
these differences increase with increasing flow driving force. To
take into account these fluctuations, we calculate the average flux,
$\bar{J}$, as
\begin{equation}
  \bar{J}=\frac{1}{n_z}\sum_{z=1}^{n_z}{J(z)}.
  \label{eq:flux}
\end{equation}
%----------------------------------
In figure \ref{f:fluxvsforce} we plot the average flux,
$\bar{J}$, as a function of applied force, $g$, in a 128$^3$ lattice
sites sample of the Bentheimer.
%----------------------------------
\begin{figure}[thb]
\centering
\includegraphics[width=0.6\textwidth]{./Fig/fluxvsforce}
\caption{
  Average flux, $\bar{J}$, as a function of applied force, $g$,
  to verify the linear behaviour predicted
  by Darcy's law (Eq.~\ref{eq:darcy}). The dotted line is the linear fit
  restricted to the data which are in the linear regime.
  Both flux and force are expressed in lattice units.
  The error bars are calculated as standard deviation over the number
  of slices along the force direction.}
\label{f:fluxvsforce}
\end{figure}


%----------------------------------
Note that the linear behaviour (with zero intercept) predicted
by Darcy's law is verified for low forcing, but deviations are
observed if the forcing is too high. %($g$=0.003 in figure \ref{f:fluxvsforce})..
For high forcing, the mean fluid velocity starts to approach the speed of sound,
$c_s$=$\frac{1}{\sqrt{3}}$,
%{\bf FIXME check this expression})
and compressibility effects become important.  The permeability can be
calculated from the slope of the plotted line in the region where
linear behaviour is observed.

%\hspace{5cm}
\begin{figure}[tb]
\centering
%\includegraphics[width=0.6\textwidth, angle=-90]{./Fig/std-collect.ps}
\includegraphics[width=0.6\textwidth, angle=-90]{./Fig/std-collect-new.ps}
\caption{Standard deviation of the average flux, $\bar{J}$, as a function of applied force, $g$.}
\label{f:stdfluxvsforce}
\end{figure}
%----------------------------------

In Figure \ref{f:stdfluxvsforce}, we show the
the standard deviation of the average flux $\bar{J}$, as a function of the force applied, $g$.
Please note that the standard deviation of the flux increases almost linearly with the force applied, for forces smaller than 0.004.
Above this value, the dependence is much stronger.
This is probably due to the following two issues: the first is the compressibility of the fluid due to the use of single relaxation time $\tau$; the second is the use of simple bounce-back boundary conditions, see Section~\ref{sec:BC}. We will investigate the relative importance of these issues in a forthcoming paper.

We find good agreement between the calculated values of the flux as a
function of applied force and the linear behaviour predicted by
Darcy's law. This is an example of how microscopic quantities, such as
the velocity field at the pore scale computed using the LB method, can be related to macroscopic, phenomenological laws.

Once the right forcing regime is found, it is sufficient to compute
the flux for one forcing level and the permeability can be directly
calculated using expression~(\ref{eq:darcy}).  For all the simulations
presented in this study we have chosen a value of the force
$g$=0.0001.

The minimum number of LB steps required to reach steady state is, in
general, system size dependent.  Using $g$=0.0001 as the driving
force, and $\tau$=1 as the BGK relaxation parameter, we have verified
that 5000 steps are sufficient for equilibration, even in the largest
sample (512$^3$) that we have considered in this work.  For this
number of steps, average quantities, such as the flux, have reached
equilibrium.  Note, however, that the local densities and velocities
fields are not yet fully equilibrated after 5000 steps, and it takes
about 8000 steps for the microscopic quantities to reach steady state.



\subsection{Effect of system size}

In this section we show how the permeability of a realistic rock
sample depends on the size of the sample, {\em i.e.} on the portion
of rock used in the simulations to calculate the flux.
Please note that the influence of the sample size was discussed by Thovert
{\it et al.} ~\cite{TYSJA01}. The influence of the boundary conditions was
discussed in this paper and also by Gerbaux {\em et al.}~\cite{GBMMVTA10}.

In this paper, the issue is
addressed by considering samples of different sizes, and by studying
the distribution of permeability as the sample size increases. The
largest size we have considered is a 512$^3$ cube, centered in the
middle of the full data set, then a 256$^3$ cube, also centered in
the middle of the full data set. The 256$^3$ sample has then been
divided into 8 and 64 cubes of size 128$^3$ and 64$^3$,
respectively. To give an idea of the physical size of these
samples, consider that a 64$^3$ cube has a linear size of
313 $\mu$m. The permeability has been computed in all these samples.
The values of permeability, $K$, as a function of void fraction
$\phi$ are reported in figure \ref{f:KvsSize}.
%-----------------
\begin{figure}[htb]
\centering
\includegraphics*[width=0.6\textwidth]{./Fig/perm_512_256_128_64.eps}
\caption{
  Permeability, $K$ (in Darcys) as a function of void fraction, $\phi$,
  for rock samples of different size. The dotted
  line connecting the points for the systems of size
  64$^3$ is only a guide to the eye.}
\label{f:KvsSize}
\end{figure}
%-----------------
Observe how the data for the smaller size we have considered (64$^3$)
are more spread-out, both in void fraction ($\phi$) and in permeability, while, as the
system size increases, the data tend to concentrate into a narrower
region. On  average, we found that the permeability increases
almost linearly with increasing void fraction, although for the
smallest samples we observed large fluctuations. 

\subsection{Effect of inlet and outlet boundary conditions}
%-------------------------------------------------------------
\begin{figure}[thb]
\centering
\subfigure[]{\includegraphics[width=0.49\textwidth]{./Fig/bc_64_2.eps}\label{f:CaseA:a}}
\subfigure[]{\includegraphics[width=0.49\textwidth]{./Fig/bc_64.eps}\label{f:CaseA:b}}
\caption{
  Permeability, K (in Darcys), as a function of void
  fraction, $\phi$, for  8 rock samples of size 64$^3$ and different
  boundary conditions (a). Figure (b) represents the same data as
  figure (a), but excluding data for the sample with largest void fraction.}
\label{f:CaseA}
\end{figure}
%---------------------------------------------------------------
The aim of this section is to investigate the influence of boundary
conditions used in LB simulations on the calculated flow field, and
on the permeability in particular.
Note that here by boundary conditions we mean the
conditions imposed on the surface of the rock sample, and not the
boundary conditions at the wall nodes, as in section \ref{sec:BC}. We
consider, for different sample sizes, the following cases:
\begin{itemize}
\item Simple periodic boundary conditions (PBC) in the three spatial directions.
\item Periodic boundary conditions with a fluid buffer 4 lattice sites wide
in the direction of the flow (direction
of the forcing, $z$-axis).
\item Mirror boundary conditions:
The system is mirrored in the direction of the flow ($z$-axis),
and then periodic boundary conditions are used.
For example, if the system size is 64$^3$ lattice sites,
then the dimensions of the mirrored system are 64x64x128 lattice sites.
\item
To study the effect of the surrounding environment and rock geometry
on the flow field, we also consider the case of a given rock sample
when taken as a sub-sample of a larger system. This means that we do
a single flow simulation for the larger system and calculate the
fluxes for the sub-samples within the larger one.
%{\bf Maddalena, is this correct? Yes it is correct}
\end{itemize}
%-----------------

Figure \ref{f:CaseA} shows the permeability for eight 64$^3$ cubes
for different boundary conditions.
As a first observation, we note that the largest changes in
permeability value with different boundary conditions are observed
for the cube with the largest void fraction and highest permeability
(last points on the right in figure \ref{f:CaseA:a}). This is the
cube where one large channel is present. If the cube is part of a
larger portion of rock, the fluid flows through different paths and
channels across the rock, and does not all flow through the large
channel. On the other hand, if this cube is taken by itself, all of the
fluid flows into the large channel, and the permeability of the
rock is higher. From figure \ref{f:CaseA:b}, where we exclude
the results for the cube discussed above, we observe that simple
PBC give lower permeability values with respect to the other
boundary conditions. This is because, by simply replicating the data
set in the direction of the flow, some channels might be
artificially closed by the replica image of the rock. The addition
of a buffer reopens these channels, hence resulting in a more
permeable medium. Also, two sets of parallel curves in $K$ versus
$\phi$ can be identified in figure \ref{f:CaseA:b},
one corresponding to the case when the cubes are isolated, with both
simple PBC and mirror PBC (open and full circles), and one for the
case when the cubes are taken as sub-samples of larger samples (open
triangles from a 256$^3$ sample and full triangles from a 512$^3$
sample).  This indicates, not surprisingly, that the flow patterns and
fluid average velocity are very much dependent on the rock geometry.
Moreover, it is also worth noticing that for all boundary conditions
considered the behavior of the permeability as a function of void
fraction is not strictly monotonous. This indicates that, besides the
void fraction, the connectivity between the pore space is also
important in determining the flow (and hence the rock permeability).

%-------------------------------------------------------------
\begin{figure}[thb]
\centering
\includegraphics[width=0.49\textwidth]{./Fig/bc_128.eps}
\caption{
  Permeability, K (in Darcys), as a function of void
  fraction, $\phi$, for  8 rock samples of size 128$^3$ and different
  boundary conditions.}
\label{f:CaseA:c}
\end{figure}
%---------------------------------------------------------------

Comparing the permeability of the 64$^3$ cubes (figure
\ref{f:CaseA:a}) with the permeability of the 128$^3$ cubes (figure
\ref{f:CaseA:c}), for different boundary conditions, it can be noted
that in the case of the larger sample the values of permeability
are less dependent on boundary conditions, and on whether the sample
is isolated or part of a larger portion of rock. This gives an
indication of the length scales over which the rock can be
considered homogeneous. Note also the narrower range of variation in void
fractions and permeability in the case of the 128$^3$ systems, as
already observed in figure \ref{f:KvsSize}.


\subsubsection{Influence of boundary conditions as a function of system size}
For this study we have considered three sub-samples of increasing
size, namely: 64$^3$, 128$^3$ and 256$^3$ voxels, respectively, all
centered at the centre of the full sample (see figure \ref{f:benth}).
%----------------------------------------------------------------------
\begin{figure}[thb]
\begin{centering}
\subfigure[256$^3$]{\includegraphics[width=0.5\textwidth]{./Fig/Layered_slice_half_y.eps}}
\subfigure[128$^3$]{\includegraphics[width=0.25\textwidth]{./Fig/128slice_half_y.eps}}
\subfigure[64$^3$]{\includegraphics[width=0.125\textwidth]{./Fig/64slice_half_y.eps}}
\caption{
  Two dimensional slices (perpendicular to the $y$-axis)
  of samples of Bentheimer sandstone.
  The slices are taken in the centre of each sample.
  Note that the larger samples contain the smaller one/ones.
  The pore space is depicted in red and the rock in black.}
\label{f:benth}
\end{centering}
\end{figure}
%----------------------------------------------------------------------
%-----------------
\begin{figure}[bth]
\centering
\includegraphics[width=0.7\textwidth]{./Fig/Size_and_BC.eps}
\caption{
  Permeability, K, in Darcys, as a function of void fraction, $\phi$,
  for rock samples of size 64$^3$, 128$^3$, 256$^3$ and
  512$^3$, for different boundary conditions.}
\label{f:Size_and_BC}
\end{figure}
%-----------------
The permeability as a function
of system size for different boundary conditions
is plotted in figure \ref{f:Size_and_BC}.
A systematic trend in the value of the permeability can be observed;
for each size, and independently of the environment, the permeability increases
in going from simple PBC, to buffer, to mirror. However, as the
sample size increases, the value of the permeability becomes less
dependent on the type of boundary conditions used. The reason for the observed trend is that, for larger samples, surface effects on the fluid flow
behaviour are smaller. For reference purposes we
also plot the permeability for the 512$^3$ sample with simple
periodic boundary conditions (star symbol).
It should also be pointed out that the data points for the smallest
system (64$^3$) show a marked dependence on the particular system
chosen, since this sample size is not large enough to capture the
length scales of the rock inhomogeneities, as also shown by the
large variations in the permeability values of figure \ref{f:KvsSize}.



\section{Displacement Propagators}
\label{propagators}

Here we report results of calculations of the
probability distribution of the solute in the
flow field computed from 3d LB simulations. The
displacement propagators from simulations can be compared with
NMR propagators. Please note that NMR has been studied by Tessier {\em et al.}~\cite{TPTA97} by means
of propagators in packed beds, both experimentally and numerically.


Following Sukop {\it et al.}~\cite{LBMODELING} and Yoshino {\it et al.} ~\cite{yoshino_2003} in their study
of dispersion porous media,
a passive scaler approach was used to solve the convecn-diffusion equation. A second distribution function $g$ with a simpler equilibrium distribution was introduced to discribe the evolution of solute concentration:

\begin{equation}
        g_i^{eq}=\omega_i \rho_s(1+3e_i \cdot u)
        \end{equation}

where the component $\rho_s$ indicate the local concentration of the solute:

\begin{equation}
        \rho_s=\sum_{i=0}^{a} g_i
        \label{rhos_m}
        \end{equation}

The following LB equation is used to compute the evolution of the distribution of solute:

\begin{equation}
        g_i(r+c_i \Delta t, t+\Delta t)-g_i(r,t)=-\frac{1}{\tau_g}[g_i(r,t)-g_i^{eq}(r,t)]
        \label{evo_g}
        \end{equation}

The diffusion coefficient $D_s$ is given by 
\begin{equation}
        D_s=\frac{1}{3}(\tau_s-0.5)
        \end{equation}

As shown in \cite{inamuro_2002} and \cite{inamuro_1997} the governing equation in macro scale can be recoverd by applying Chapman-Enskog expansion to equation \ref{evo_g} and \ref{rhos_m}:


\begin{equation}
        \frac{\partial \rho_s}{\partial t}=-\nu \nabla \rho_s + D_s \Delta \rho_s
        \end{equation}


This approach is only able to calculate the evolution of the solute concentration, to compute the tracer particles evolution, we initialize the solute concentration with a peak distribution on lattice $i$ in X direction:

\begin{equation}
        \rho_s(r)=1 \quad r_x=i; \quad \rho_s(r)=0 \quad \text{ elsewise}
        \label{ini_rhos}
        \end{equation}

With this initial condition, we can simulate the evolution of tracer particles from location $r_x=i$. If uniformed ditribute particles are to computed, we can add more "solutes" to calculate the displacement of tracer particles from varying starting point. This will increase the computing time and memory use significantly; however, if the code has been optimized well, the computing time is still acceptable.  





\subsection{Results}

In Figure~\ref{f:zeta-t}, we show the mean displacement $\left<
\zeta \right>$ normalised by the nominal mean displacement $\left<
\zeta \right>_0$ = $\bar{v} T$ as a function of simulation time
$T$ for a system of size 128$^3$.

%----------------------------------
\begin{figure}[tb]
\centering
\includegraphics[width=0.5\textwidth]{./Fig/ratio_mean_disp_nominal_mean_disp.eps}
\caption{
  Mean displacement $\left< \zeta \right>$ normalised by the
  nominal mean displacement $\left< \zeta \right>_0$ as a function of
  evolution time $T$. $\bar{v}=0.4 mm/s$}
\label{f:zeta-t}
\end{figure}
%----------------------------------

$\left< \zeta \right>_0$ is the mean displacement along the mean
flow direction of a fluid flowing with the constant mean velocity $\bar{v}$
during the time evolution interval $T$. We observe
that the value of $\zeta / \left< \zeta \right>_0$ converge to 0.996 which is very close to 1 as we expected for long time hehaviour. The small deviation of 0.004 is prbably due to 1) numerical
rounding off errors in the concentration calculation.

In Figure~\ref{f:Pz-zeta}, we show the probability distribution of
displacements $P\left(\zeta \right)$ as a function of the
displacement $\zeta$.
%----------------------------------
\begin{figure}[htb]
\centering
\includegraphics[width=0.8\textwidth]{./Fig/Bentheimer_512_concentration_displacement_plot.eps}
\caption{
  Probability distribution of displacements $P\left(\zeta
   \right)$ as a function of the displacement $\zeta$ for different simulation
    times. $\bar{v}=0.4 mm/s$.}
\label{f:Pz-zeta}
\end{figure}
%----------------------------------
From this figure, we observe that the displacement propagator
distribution evolves in time. The propagator distributions show dispersive behaviour, in
qualitative agreement with recent experiments by Scheven {\em et
al.}~\cite{bib:schevenpf2006}.

%In Figure~\ref{f:Pz-zeta-norm}, we show the time evolution of the
%normalised probability distribution of displacements
%$\left< \zeta \right>_0 P\left(\zeta \right)$ as a function of normalised
%displacement $\zeta / \left< \zeta \right>_0$.
%----------------------------------
%\begin{figure}[htb]
%\centering
%\includegraphics[angle=270,width=0.5\textwidth]{./Fig/propagator/propagators_normalised-l.eps}
%\caption{
%  Normalised probability distribution of displacements
%  $ \left< \zeta \right>_0 P\left(\zeta \right) $ as a function of
%  normalised displacement $\zeta / \left< \zeta \right>_0$.}
%\label{f:Pz-zeta-norm}
%\end{figure}
%----------------------------------
%We observe that the distribution shifts to higher values of $\zeta /
%\left< \zeta \right>_0$ with increasing evolution time, in
%qualitative agreement with the experiments by Scheven {\em et
%al.}~\cite{bib:schevenpf2006} for a Bentheimer sandstone.


%%%%%%% Two-phase flow %%%%%%%
%______Parameter study for binary fluid
\section{Phase Separation of binary mixtures}

In this section we present results of immiscible fluid flow
using the Shan-Chen model~\cite{bib:shan-chen} discussed in
section~\ref{Sec:LatticeBoltzmann}. We first study phase
separation and wettability effects. Then we describe
calculations of relative permeabilities in the Bentheimer rock sample.

The phase separation of a binary immiscible mixture in LB simulations depends on a
number of parameters, particularly  the value of the
coupling constant  $g_{\sigma \bar{\sigma}}$ between the two phases
(see equation \ref{eq:colour-colour}). Previous studies (see e.g.~\cite{bib:chin-boek-coveney})
have shown that the surface tension increases linearly with increasing coupling constant.
LB simulations of critical spinodal decomposition in
binary immiscible fluids~\cite{bib:gonzalez-nekovee-coveney} were also
carried out to determine the properties of a binary mixture as
function of different parameter sets.
Different wettabilities are implemented in our model
by assigning a particular density of the wetting phase to the solid sites (the surface).
This density is not evolving with the LB equation,
but exerts a force on the neighboring fluid sites
of the opposite phase, with the same
coupling constant as the fluid-fluid interaction.
For example, for water wettability, if the bulk density
of the water phase is $\rho=0.8$, then a density between 0 (non wet)
and 0.8 (maximal wettability) can be assigned to the solid sites.
The value of the wettability is given by the solid site density divided by the
wetting fluid bulk density.

To verify that our implementation of two-phase flow simulations with the LB method is able to
model a wide range of wettability situations, we used a simple
model of a box containing two immiscible fluids above a surface with
variable wettability.
The non-wetting fluid
is initially placed in a cubic configuration at the center of the box and in contact
with the surface, then a LB simulation is performed until
the non-wetting phase has reached its equilibrium shape.
The results are shown in figure~\ref{fig:dropletwet} for oil
concentration of 25$\%$ (the oil phase is colored in red).
When the surface is non wet (wet=0), the non-wetting (oil) phase spreads
over the entire surface.
As the water wettability of the surface is increased, the non-wetting phase starts to detach from the surface,
and forms a contact angle with it. This contact angle
can be defined as the angle between the wet surface and a plane tangent to the non-wetting
phase droplet where the droplet starts to detach from the surface,
as illustrated in figure~\ref{fig:contactangle}.
This contact angle increases with increasing wettability.
At high wettability ($\ge$0.6)  the
non-wetting phase is completely detached from the
surface and forms a spherical bubble.

\begin{figure}[htbp]
\begin{center}
\includegraphics[width=0.2\textwidth]{./Fig/contactangle.eps}
\end{center}
\caption{
  Schematic illustration of how the contact angle $\theta$ is defined.
}
\label{fig:contactangle}
\end{figure}
%--------------

\begin{figure}[htbp]
\begin{center}
\subfigure[wet=0.0]{\includegraphics[width=0.2\textwidth,angle=180]{./Fig/PipeFr-19-10/snap_ps_rho0.8_pipe_wet0.0.eps}}
\subfigure[wet=0.1]{\includegraphics[width=0.2\textwidth,angle=180]{./Fig/PipeFr-19-10/snap_ps_rho0.8_pipe_wet0.1.eps}}
\subfigure[wet=0.2]{\includegraphics[width=0.2\textwidth,angle=180]{./Fig/PipeFr-19-10/snap_ps_rho0.8_pipe_wet0.2.eps}}
\subfigure[wet=0.3]{\includegraphics[width=0.2\textwidth,angle=180]{./Fig/PipeFr-19-10/snap_ps_rho0.8_pipe_wet0.3.eps}}
\subfigure[wet=0.4]{\includegraphics[width=0.2\textwidth,angle=180]{./Fig/PipeFr-19-10/snap_ps_rho0.8_pipe_wet0.4.eps}}
\subfigure[wet=0.5]{\includegraphics[width=0.2\textwidth,angle=180]{./Fig/PipeFr-19-10/snap_ps_rho0.8_pipe_wet0.5.eps}}
\subfigure[wet=0.6]{\includegraphics[width=0.2\textwidth,angle=180]{./Fig/PipeFr-19-10/snap_ps_rho0.8_pipe_wet0.6.eps}}
\subfigure[wet=0.8]{\includegraphics[width=0.2\textwidth,angle=180]{./Fig/PipeFr-19-10/snap_ps_rho0.8_pipe_wet0.8.eps}}
\end{center}
\caption{
  Equilibrium configurations (as 2d slices taken at the center of the 3d system)
  of a binary immiscible fluid mixture
  as a function of increasing water wettability of the bottom surface.
  Wetting phase (water) is depicted in blue and non-wetting phase (oil) in red.
  Note the finite size effects in panel (c) where
  the detaching droplet of oil touches the lateral boundaries of the box.
  In this case it is not possible
  to determine the contact angle.
  However, there are no finite size effects for all the other
  values of wettability.
}
\label{fig:dropletwet}
\end{figure}
%--------------
%--------------



%-----------------------------
%\begin{table}[htbp] \centering
%\caption{Binary mixtures phase separation.  $\psi$=$\rho$, bdist=0, g$_{\rm br}$=0.06}
%\label{t:1}
%\begin{tabular}{c|cccc|l}
%              &            &             &           &        &   \\
%              &   0.6      &      0.7    &    0.8    &    0.9 &   \\
%\hline
%       &    No      &     No      &   No      &   No   &  bulk  \\
% 0.1   &    No      &     No      &   Partial &   Yes  &  rock (wet=0)  \\
%       &    Partial &    Partial  &   Partial &   Yes  &  rock (wet=-0.4)  \\
%\hline
%       &    No      &    Yes      &   Yes     &  nan   & bulk    \\
% 0.2   &    No      &    Partial  &   Yes     &  nan   & rock (wet=0)   \\
%       &    Partial &    Yes      &   Yes     &  nan   & rock (wet=-0.4)   \\
%\hline
%       &    Yes     &    Yes      &   Yes     &  nan   & bulk    \\
% 0.3   &    Partial &    Yes      &   Yes     &  nan   & rock (wet=0)    \\
%       &    Yes     &    Yes      &   Yes     &  nan   & rock (wet=-0.4)    \\
%\hline
%       &    Yes     &    Yes      &   Yes     &  nan   & bulk    \\
% 0.4   &    Yes     &    Yes      &   Yes     &  nan   & rock (wet=0)    \\
%       &    Yes     &    Yes      &   Yes     &  nan   & rock (wet=-0.4)     \\
%\hline
%       &    Yes     &    Yes      &   Yes     &   nan  & bulk   \\
% 0.5   &    Yes     &    Yes      &   Yes     &   nan  & rock (wet=0)   \\
%       &    Yes     &    Yes      &   nan     &   nan  & rock (wet=-0.4)    \\
%\hline
%       &  --------  & --------    &  -------- &  -------- &   \\
% 0.6   &  --------  & --------    &  -------- &  -------- &  \\
%       &    Yes     &    Yes      &  nan      &    nan    & rock (wet=-0.4)  \\
%\hline
%       & --------   & --------    & --------  &  -------- &   \\
% 0.7   & --------   & --------    & --------  &  -------- &  \\
%       &    Yes     &    Yes      &   nan     &    nan    & rock (wet=-0.4)  \\
%\hline
%       & --------   &  --------   & --------  &  -------- &   \\
% 0.8   & --------   &  --------   & --------  &  -------- &  \\
%       &    No      &    Yes      &   nan     &     nan   &  rock (wet=-0.4) \\
%\hline
%       & --------   &  --------   &  -------- &  -------- &   \\
% 0.9   & --------   &  --------   &  -------- &  -------- &  \\
%       &    No      &    No       &    No     &     No    & rock (wet=-0.4)  \\
%\hline
%\end{tabular}
%\end{table}




%\begin{table}[htbp] \centering
%\caption{Binary mixtures phase separation. $\psi$=$\rho$, bdist=2, $\rho$=0.8, g$_{\rm br}$=0.06}
%\label{t:2}
%\begin{tabular}{c|cc|l}
%        &                      &                          &             \\
%        &   g$_{\rm accn}$=0   &   g$_{\rm accn}$=0.0005  &             \\
%\hline
% 0.1    &  Yes (not complete)  &            ?             & rock (wet=-0.4)          \\
%\hline
% 0.2    &      Yes             &        Yes               & rock (wet=-0.4)          \\
%\hline
% 0.3    &      Yes             &        Yes               & rock (wet=-0.4)          \\
%\hline
% 0.4    &     Yes              &        Yes               & rock (wet=-0.4)          \\
%\hline
% 0.5    &     Yes              &        Yes               & rock (wet=-0.4)          \\
%\hline
% 0.6    &     Yes              &        Yes               & rock (wet=-0.4)          \\
%\hline
% 0.7    &     Yes              &        Yes               & rock (wet=-0.4)          \\
%\hline
% 0.8    &     Yes              &        Yes               & rock (wet=-0.4)          \\
%\hline
% 0.9    &     No               &        No                & rock (wet=-0.4)          \\
%\hline
%\end{tabular}
%\end{table}



%______Relative permeabilities
\section{Relative permeabilities}
\label{relperm}

We carried out an extensive parameter study to find the
value of the total density that gives both phase separation and
stability. Using the optimal value of 0.8, we calculated the relative
permeabilities for both oil and water, applying a driving force
to either oil, or water, or both.
The initial condition for the distribution of fluids was a
homogeneous saturation of the rock, where the saturation at each
node is fixed. This means that, at each lattice site, the fractions
of oil and water are fixed. In all cases, we run the simulations
until a steady state is achieved. At steady state the two fluids
are phase separated.
The results are shown in Figure~\ref{f:relperm}. In all
cases, the rock is considered to be water-wet. We observe that the
results for the non-wetting fluid (oil) are independent of the
details of the forcing. For the case of the wetting fluid (water),
the results depend on whether we force water or both water and oil.
In a future study, we will present a more detailed analysis including comparison with different models for immiscible fluids~\cite{bib:keehm-2004}.

%----------------------------------
\begin{figure}[htb]
\centering
\includegraphics[width=0.6\textwidth,angle=-90]{./Fig/relperm_all_cases.ps}
\caption{
  Relative permeabilities for wetting (water) and non-wetting (oil) phase
  in a sample of Bentheimer sandstone of 128$^3$ voxels.
}
\label{f:relperm}
\end{figure}
%----------------------------------


%______Conclusions
\section{Conclusions}
We have used lattice-Boltzmann simulations to calculate the
single-phase permeability of a realistic three-dimensional sample
of Bentheimer sandstone. By systematically studying
the effect of system size and periodic boundary conditions, we
were able to validate Darcy's law from the linear dependence of the flux on the
body force exterted. Our results show that
the values of the permeability
measurements as a function of porosity tend to concentrate in a
narrower region of the porosity as the system size of the simulated
rock sample increases.
By investigating the effect of different boundary conditions on rock
samples of different size, we were able to infer the length scales
over which the rock can be considered homogenous.
Calculations of displacement propagators using the LB
flow field show dispersive behaviour, in agreement with NMR
experiments. Finally, we have presented some results on
relative permeabilities in Bentheimer sandstone for
a binary immiscible fluid at different forcing conditions.
The latter point is of great interest for practical applications, and
deserves further studies.

\section*{Acknowledgments}

We are grateful to Matt Harvey  for help with the propagator calculations and
to John Crawshaw for helpful discussions.
The research was supported by
the EPSRC RealityGrid grants GR/R67699/01 and GR/R67699/02.
To carry out the most computationally expensive simulations presented in this paper
resources on a Cray XT3 (``LeMieux'')
at Pittsburgh PSC HPC facility (USA) were used.
% \input{acknowledgments}

%______Bibliography
\bibliographystyle{unsrt}
%\bibliographystyle{plain} % alphabetical
\bibliography{refs}

\end{document}
